/*
 * Academic License - for use in teaching, academic research, and meeting
 * course requirements at degree granting institutions only.  Not for
 * government, commercial, or other organizational use.
 *
 * Modalmass_estimation_newchangefinal1219.c
 *
 * Code generation for function 'Modalmass_estimation_newchangefinal1219'
 *
 */

/* Include files */
#include "Modalmass_estimation_newchangefinal1219.h"
#include "emxutil.h"
#include "Modalmass_estimation_newchangefinal1219_rtwutil.h"
#include "types.h"
#include "fsolve.h"
#include "minOrMax.h"
#include "ppval.h"
#include "rms.h"
#include "rt_nonfinite.h"
#include "sparse.h"
#include "sparse1.h"
#include "spline.h"
#include "wavedec.h"
#include <math.h>
#include <string.h>

/* Function Declarations */
static void binary_expand_op(emxArray_real_T *in1, const emxArray_int8_T *in2,
                             const emxArray_int8_T *in3,
                             const emxArray_real_T *in4,
                             const emxArray_int32_T *in5,
                             const emxArray_int32_T *in6, int in7, int in8);

/* Function Definitions */
static void binary_expand_op(emxArray_real_T *in1, const emxArray_int8_T *in2,
                             const emxArray_int8_T *in3,
                             const emxArray_real_T *in4,
                             const emxArray_int32_T *in5,
                             const emxArray_int32_T *in6, int in7, int in8)
{
  emxArray_real_T *b_in1;
  double *b_in1_data;
  double *in1_data;
  int aux_0_1;
  int aux_1_1;
  int aux_2_1;
  int b_loop_ub;
  int i;
  int i1;
  int loop_ub;
  int stride_0_0;
  int stride_0_1;
  int stride_1_0;
  int stride_1_1;
  int stride_2_0;
  int stride_2_1;
  const signed char *in2_data;
  const signed char *in3_data;
  in3_data = in3->data;
  in2_data = in2->data;
  in1_data = in1->data;
  emxInit_real_T(&b_in1, 2);
  if (in3->size[0] == 1) {
    if (in2->size[0] == 1) {
      loop_ub = in1->size[0];
    } else {
      loop_ub = in2->size[0];
    }
  } else {
    loop_ub = in3->size[0];
  }
  i = b_in1->size[0] * b_in1->size[1];
  b_in1->size[0] = loop_ub;
  if (in3->size[1] == 1) {
    if (in2->size[1] == 1) {
      b_loop_ub = in1->size[1];
    } else {
      b_loop_ub = in2->size[1];
    }
  } else {
    b_loop_ub = in3->size[1];
  }
  b_in1->size[1] = b_loop_ub;
  emxEnsureCapacity_real_T(b_in1, i);
  b_in1_data = b_in1->data;
  stride_0_0 = (in1->size[0] != 1);
  stride_0_1 = (in1->size[1] != 1);
  stride_1_0 = (in2->size[0] != 1);
  stride_1_1 = (in2->size[1] != 1);
  stride_2_0 = (in3->size[0] != 1);
  stride_2_1 = (in3->size[1] != 1);
  aux_0_1 = 0;
  aux_1_1 = 0;
  aux_2_1 = 0;
  for (i = 0; i < b_loop_ub; i++) {
    for (i1 = 0; i1 < loop_ub; i1++) {
      b_in1_data[i1 + b_in1->size[0] * i] =
          (in1_data[i1 * stride_0_0 + in1->size[0] * aux_0_1] +
           (double)in2_data[i1 * stride_1_0 + in2->size[0] * aux_1_1]) +
          (double)in3_data[i1 * stride_2_0 + in3->size[0] * aux_2_1];
    }
    aux_2_1 += stride_2_1;
    aux_1_1 += stride_1_1;
    aux_0_1 += stride_0_1;
  }
  sparse_plus(b_in1, in4, in5, in6, in7, in8, in1);
  emxFree_real_T(&b_in1);
}

void Modalmass_estimation_newchangefinal1219(
    const double person[3679], double structure[7298], const double f[3],
    double weight, double t0, double d, double n, double L, double B,
    char stru_system, double original_Modalmass[3], double Modalmass[3],
    double *original_eachmass, double *eachmass, double *original_k, double *k)
{
  static const double a[101] = {0.0,
                                0.01,
                                0.02,
                                0.03,
                                0.04,
                                0.05,
                                0.06,
                                0.07,
                                0.08,
                                0.09,
                                0.1,
                                0.11,
                                0.12,
                                0.13,
                                0.14,
                                0.15,
                                0.16,
                                0.17,
                                0.18,
                                0.19,
                                0.2,
                                0.21,
                                0.22,
                                0.23,
                                0.24,
                                0.25,
                                0.26,
                                0.27,
                                0.28,
                                0.29,
                                0.3,
                                0.31,
                                0.32,
                                0.33,
                                0.34,
                                0.35000000000000003,
                                0.36,
                                0.37,
                                0.38,
                                0.39,
                                0.4,
                                0.41000000000000003,
                                0.42,
                                0.43,
                                0.44,
                                0.45,
                                0.46,
                                0.47000000000000003,
                                0.48,
                                0.49,
                                0.5,
                                0.51,
                                0.52,
                                0.53,
                                0.54,
                                0.55,
                                0.56,
                                0.57000000000000006,
                                0.58000000000000007,
                                0.59,
                                0.6,
                                0.61,
                                0.62,
                                0.63,
                                0.64,
                                0.64999999999999991,
                                0.65999999999999992,
                                0.66999999999999993,
                                0.67999999999999994,
                                0.69,
                                0.7,
                                0.71,
                                0.72,
                                0.73,
                                0.74,
                                0.75,
                                0.76,
                                0.77,
                                0.78,
                                0.79,
                                0.8,
                                0.81,
                                0.82000000000000006,
                                0.83,
                                0.84,
                                0.85,
                                0.86,
                                0.87,
                                0.88,
                                0.89,
                                0.9,
                                0.91,
                                0.92,
                                0.92999999999999994,
                                0.94,
                                0.95,
                                0.96,
                                0.97,
                                0.98,
                                0.99,
                                1.0};
  static double b_value[10000];
  emxArray_int32_T *t2_colidx;
  emxArray_int32_T *t2_rowidx;
  emxArray_int8_T *c_d;
  emxArray_int8_T *d_d;
  emxArray_real_T *b_C;
  emxArray_real_T *b_I;
  emxArray_real_T *c;
  emxArray_real_T *c_I;
  emxArray_real_T *t2_d;
  emxArray_real_T *y_val;
  emxArray_real_T *y_val2;
  emxArray_real_T *y_val3;
  double A[2000];
  double F[2000];
  double a__3[2000];
  double acol_data[2000];
  double b_acol_data[2000];
  double t1_coefs[400];
  double phi[303];
  double b_y[101];
  double c_y[101];
  double y[101];
  double l[7];
  double x[4];
  double C;
  double Di;
  double Ki;
  double a_tmp;
  double alpha_0;
  double alpha_idx_0;
  double alpha_idx_1;
  double alpha_idx_2;
  double b_d;
  double b_expl_temp;
  double bsum;
  double original_Modal_Mass1;
  double original_mass;
  double *C_data;
  double *I_data;
  double *c_data;
  double *y_val2_data;
  double *y_val3_data;
  double *y_val_data;
  int acol_size;
  int b_k;
  int b_loop_ub_tmp;
  int c_i;
  int c_loop_ub_tmp;
  int i;
  int loop_ub_tmp;
  int m0;
  int nv_tmp_tmp_tmp;
  short iv[2000];
  signed char *b_d_data;
  signed char *d_data;
  boolean_T b;
  /*  输入参数: */
  /*  - person: 人体侧加速度(g) */
  /*  - structure: 结构侧加速度(g) */
  /*  - f: 基频(Hz) */
  /*  - weight: 重量 (kg) */
  /*  - d: 阻尼比 */
  /*  - n: 层数 (integer) */
  /*  - L: 结构长度 (m) */
  /*  - B: 结构宽度 (m) */
  /*  - stru_system: 结构类型 ('框架结构', '剪力墙结构', '框架剪力墙结构') */
  /*  参数设置 */
  bsum = t0 * 100.0 + 1.0;
  Ki = weight * 9.81 * (0.23 * f[0] + 0.339);
  for (i = 0; i < 2000; i++) {
    short b_i;
    b_i = (short)(bsum + (double)i);
    iv[i] = b_i;
    F[i] = Ki * person[b_i - 1];
  }
  Ki = F[0];
  for (b_k = 0; b_k < 1023; b_k++) {
    Ki += F[b_k + 1];
  }
  bsum = F[1024];
  for (b_k = 0; b_k < 975; b_k++) {
    bsum += F[b_k + 1025];
  }
  Ki += bsum;
  Ki /= 2000.0;
  for (i = 0; i < 2000; i++) {
    F[i] -= Ki;
  }
  for (i = 0; i < 7298; i++) {
    structure[i] *= 9.81;
  }
  Ki = structure[iv[0] - 1];
  for (b_k = 0; b_k < 1023; b_k++) {
    Ki += structure[iv[b_k + 1] - 1];
  }
  bsum = structure[iv[1024] - 1];
  for (b_k = 0; b_k < 975; b_k++) {
    bsum += structure[iv[b_k + 1025] - 1];
  }
  Ki += bsum;
  Ki /= 2000.0;
  /*  小波变换参数 */
  for (i = 0; i < 2000; i++) {
    acol_data[i] = structure[iv[i] - 1] - Ki;
  }
  emxInit_real_T(&c, 1);
  wavedec(acol_data, c, l);
  c_data = c->data;
  acol_size = (int)l[0];
  for (b_k = 0; b_k < acol_size; b_k++) {
    b_acol_data[b_k] = c_data[b_k];
  }
  /*  初始化向量 */
  a_tmp = 6.2831853071795862 * f[0];
  for (i = 0; i < 10000; i++) {
    m0 = 1000 * (i + 1);
    /*  输入 */
    bsum = (double)m0 * (a_tmp * a_tmp);
    C = 2.0 * d * (6.2831853071795862 * f[0]) * (double)m0;
    /*  使用 Newmark_ode 求解运动方程 */
    /*  Newmark方法求解ODE */
    memset(&acol_data[0], 0, 2000U * sizeof(double));
    memset(&a__3[0], 0, 2000U * sizeof(double));
    memset(&A[0], 0, 2000U * sizeof(double));
    A[0] = ((F[0] - bsum * 0.0) - C * 0.0) / (double)m0;
    Ki = (60000.0 * (double)m0 + 300.0 * C) + bsum;
    acol_data[0] = 0.0;
    a__3[0] = 0.0;
    for (c_i = 0; c_i < 1999; c_i++) {
      b_d = acol_data[c_i];
      b_expl_temp = a__3[c_i];
      original_mass = A[c_i];
      Di = ((F[c_i + 1] + (double)m0 * ((60000.0 * b_d + 600.0 * b_expl_temp) +
                                        2.0 * original_mass)) +
            C * ((300.0 * b_d + 2.0 * b_expl_temp) + 0.005 * original_mass)) /
           Ki;
      bsum = (60000.0 * (Di - b_d) - 600.0 * b_expl_temp) - 2.0 * original_mass;
      a__3[c_i + 1] = (a__3[c_i] + 0.005 * original_mass) + 0.005 * bsum;
      acol_data[c_i + 1] = Di;
      A[c_i + 1] = bsum;
    }
    /*  小波变换 */
    wavedec(A, c, l);
    c_data = c->data;
    m0 = (int)l[0];
    for (b_k = 0; b_k < m0; b_k++) {
      acol_data[b_k] = c_data[b_k];
    }
    /*  计算差值乘以权重 */
    b_value[i] = fabs(rms(acol_data, (int)l[0]) - rms(b_acol_data, acol_size));
  }
  double a__5[4];
  char expl_temp[19];
  minimum(b_value, &m0);
  original_Modal_Mass1 = 1000.0 * (double)m0;
  /*  使用 root3d 函数求解 */
  /*  调用 fsolve */
  fsolve(f, x, a__5, &Ki, &Di, &b_expl_temp, &original_mass, expl_temp);
  /*  计算 Modal Mass */
  emxInit_real_T(&y_val, 2);
  i = y_val->size[0] * y_val->size[1];
  y_val->size[0] = 1;
  loop_ub_tmp = (int)n;
  y_val->size[1] = (int)n;
  emxEnsureCapacity_real_T(y_val, i);
  y_val_data = y_val->data;
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val_data[i] = 0.0;
  }
  alpha_0 = x[3];
  /*  子函数 calculate_phi */
  for (i = 0; i < 3; i++) {
    double b_etai_tmp;
    double b_x;
    double etai;
    double etai_tmp;
    b_x = x[i];
    etai_tmp = sin(b_x);
    Ki = cos(b_x);
    bsum = b_x * b_x;
    b_expl_temp = alpha_0 * alpha_0 + bsum;
    b_etai_tmp = sqrt(b_expl_temp);
    Di = sinh(b_etai_tmp);
    original_mass = cosh(b_etai_tmp);
    etai = (bsum * etai_tmp + b_x * b_etai_tmp * Di) /
           (bsum * Ki + b_expl_temp * original_mass);
    C = b_x * rt_powd_snf(b_expl_temp, -0.5);
    etai_tmp = (etai_tmp - C * Di) + etai * (original_mass - Ki);
    for (b_k = 0; b_k < 101; b_k++) {
      b_d = a[b_k];
      b_expl_temp = b_d * b_etai_tmp;
      original_mass = sinh(b_expl_temp);
      y[b_k] = original_mass;
      b_expl_temp = cosh(b_expl_temp);
      b_y[b_k] = b_expl_temp;
      bsum = b_x * b_d;
      phi[i + 3 * b_k] =
          ((sin(bsum) - C * original_mass) + etai * (b_expl_temp - cos(bsum))) /
          etai_tmp;
    }
  }
  /*  计算质量矩阵 */
  if ((int)n - 1 >= 0) {
    for (i = 0; i < 101; i++) {
      y[i] = a[i] * n;
      b_y[i] = phi[3 * i];
    }
  }
  if (loop_ub_tmp - 1 >= 0) {
    spline(y, b_y, c_y, t1_coefs);
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val_data[i] = ppval(c_y, t1_coefs, (double)i + 1.0);
  }
  i = c->size[0];
  c->size[0] = y_val->size[1];
  emxEnsureCapacity_real_T(c, i);
  c_data = c->data;
  m0 = y_val->size[1];
  for (i = 0; i < m0; i++) {
    c_data[i] = y_val_data[i];
  }
  emxInit_real_T(&b_I, 2);
  i = b_I->size[0] * b_I->size[1];
  b_I->size[0] = (int)n;
  b_I->size[1] = (int)n;
  emxEnsureCapacity_real_T(b_I, i);
  I_data = b_I->data;
  b_loop_ub_tmp = (int)n * (int)n;
  for (i = 0; i < b_loop_ub_tmp; i++) {
    I_data[i] = 0.0;
  }
  if ((int)n > 0) {
    for (b_k = 0; b_k < loop_ub_tmp; b_k++) {
      I_data[b_k + b_I->size[0] * b_k] = 1.0;
    }
  }
  acol_size = y_val->size[1];
  c_i = b_I->size[1];
  emxInit_real_T(&b_C, 2);
  i = b_C->size[0] * b_C->size[1];
  b_C->size[0] = 1;
  b_C->size[1] = b_I->size[1];
  emxEnsureCapacity_real_T(b_C, i);
  C_data = b_C->data;
  for (i = 0; i < c_i; i++) {
    m0 = i * b_I->size[0];
    C_data[i] = 0.0;
    for (b_k = 0; b_k < acol_size; b_k++) {
      C_data[i] += y_val_data[b_k] * I_data[m0 + b_k];
    }
  }
  C = 0.0;
  m0 = b_C->size[1];
  for (i = 0; i < m0; i++) {
    C += C_data[i] * c_data[i];
  }
  original_mass = original_Modal_Mass1 * n / C;
  if (stru_system == 'A') {
    m0 = 0;
  } else if (stru_system == 'B') {
    m0 = 1;
  } else if (stru_system == 'C') {
    m0 = 2;
  } else {
    m0 = -1;
  }
  switch (m0) {
  case 0:
    alpha_idx_0 = 1.02;
    alpha_idx_1 = 0.94;
    alpha_idx_2 = 1.31;
    break;
  case 1:
    alpha_idx_0 = 1.09;
    alpha_idx_1 = 0.84;
    alpha_idx_2 = 1.52;
    break;
  case 2:
    alpha_idx_0 = 0.85;
    alpha_idx_1 = 0.73;
    alpha_idx_2 = 2.29;
    break;
  }
  *original_eachmass = original_mass / n;
  emxInit_real_T(&y_val2, 2);
  i = y_val2->size[0] * y_val2->size[1];
  y_val2->size[0] = 1;
  y_val2->size[1] = (int)n;
  emxEnsureCapacity_real_T(y_val2, i);
  y_val2_data = y_val2->data;
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val2_data[i] = 0.0;
  }
  if ((int)n - 1 >= 0) {
    for (i = 0; i < 101; i++) {
      y[i] = a[i] * n;
      b_y[i] = phi[3 * i + 1];
    }
  }
  if (loop_ub_tmp - 1 >= 0) {
    spline(y, b_y, c_y, t1_coefs);
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val2_data[i] = ppval(c_y, t1_coefs, (double)i + 1.0);
  }
  emxInit_real_T(&y_val3, 2);
  i = y_val3->size[0] * y_val3->size[1];
  y_val3->size[0] = 1;
  y_val3->size[1] = (int)n;
  emxEnsureCapacity_real_T(y_val3, i);
  y_val3_data = y_val3->data;
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val3_data[i] = 0.0;
  }
  if ((int)n - 1 >= 0) {
    for (i = 0; i < 101; i++) {
      y[i] = a[i] * n;
      b_y[i] = phi[3 * i + 2];
    }
  }
  if (loop_ub_tmp - 1 >= 0) {
    spline(y, b_y, c_y, t1_coefs);
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val3_data[i] = ppval(c_y, t1_coefs, (double)i + 1.0);
  }
  /*  计算刚度矩阵和 Modal Mass */
  i = b_I->size[0] * b_I->size[1];
  b_I->size[0] = (int)n;
  b_I->size[1] = (int)n;
  emxEnsureCapacity_real_T(b_I, i);
  I_data = b_I->data;
  for (i = 0; i < b_loop_ub_tmp; i++) {
    I_data[i] = 0.0;
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    I_data[i + b_I->size[0] * i] = 2.0;
  }
  nv_tmp_tmp_tmp = (int)(n - 1.0);
  emxInit_int8_T(&c_d);
  i = c_d->size[0] * c_d->size[1];
  c_d->size[0] = (int)(n - 1.0) + 1;
  c_d->size[1] = (int)(n - 1.0) + 1;
  emxEnsureCapacity_int8_T(c_d, i);
  d_data = c_d->data;
  c_loop_ub_tmp = ((int)(n - 1.0) + 1) * ((int)(n - 1.0) + 1);
  for (i = 0; i < c_loop_ub_tmp; i++) {
    d_data[i] = 0;
  }
  for (i = 0; i < nv_tmp_tmp_tmp; i++) {
    d_data[(i + c_d->size[0] * i) + 1] = -1;
  }
  emxInit_int8_T(&d_d);
  i = d_d->size[0] * d_d->size[1];
  d_d->size[0] = (int)(n - 1.0) + 1;
  d_d->size[1] = (int)(n - 1.0) + 1;
  emxEnsureCapacity_int8_T(d_d, i);
  b_d_data = d_d->data;
  for (i = 0; i < c_loop_ub_tmp; i++) {
    b_d_data[i] = 0;
  }
  for (i = 0; i < nv_tmp_tmp_tmp; i++) {
    b_d_data[i + d_d->size[0] * (i + 1)] = -1;
  }
  emxInit_real_T(&t2_d, 1);
  emxInit_int32_T(&t2_colidx);
  emxInit_int32_T(&t2_rowidx);
  m0 = sparse(n, n, n, n, t2_d, t2_colidx, t2_rowidx, &acol_size);
  if (b_I->size[0] == 1) {
    i = c_d->size[0];
  } else {
    i = b_I->size[0];
  }
  b = ((b_I->size[0] == c_d->size[0]) && (b_I->size[1] == c_d->size[1]) &&
       (i == d_d->size[0]) && (i == d_d->size[1]));
  emxInit_real_T(&c_I, 2);
  if (b) {
    i = c_I->size[0] * c_I->size[1];
    c_I->size[0] = b_I->size[0];
    c_I->size[1] = b_I->size[1];
    emxEnsureCapacity_real_T(c_I, i);
    C_data = c_I->data;
    for (i = 0; i < b_loop_ub_tmp; i++) {
      C_data[i] = (I_data[i] + (double)d_data[i]) + (double)b_d_data[i];
    }
    sparse_plus(c_I, t2_d, t2_colidx, t2_rowidx, m0, acol_size, b_I);
    I_data = b_I->data;
  } else {
    binary_expand_op(b_I, c_d, d_d, t2_d, t2_colidx, t2_rowidx, m0, acol_size);
    I_data = b_I->data;
  }
  acol_size = y_val->size[1];
  c_i = b_I->size[1];
  i = b_C->size[0] * b_C->size[1];
  b_C->size[0] = 1;
  b_C->size[1] = b_I->size[1];
  emxEnsureCapacity_real_T(b_C, i);
  C_data = b_C->data;
  for (i = 0; i < c_i; i++) {
    m0 = i * b_I->size[0];
    C_data[i] = 0.0;
    for (b_k = 0; b_k < acol_size; b_k++) {
      C_data[i] += y_val_data[b_k] * I_data[m0 + b_k];
    }
  }
  bsum = 0.0;
  m0 = b_C->size[1];
  for (i = 0; i < m0; i++) {
    bsum += C_data[i] * c_data[i];
  }
  Ki = a_tmp * a_tmp;
  *original_k = Ki * original_Modal_Mass1 / bsum;
  /*  返回 Modal Mass 和矩阵 */
  bsum = 0.0;
  m0 = y_val2->size[1];
  for (i = 0; i < m0; i++) {
    b_d = y_val2_data[i];
    bsum += b_d * *original_eachmass * b_d;
  }
  b_expl_temp = 0.0;
  m0 = y_val3->size[1];
  for (i = 0; i < m0; i++) {
    b_d = y_val3_data[i];
    b_expl_temp += b_d * *original_eachmass * b_d;
  }
  original_Modalmass[0] = original_Modal_Mass1;
  original_Modalmass[1] = bsum;
  original_Modalmass[2] = b_expl_temp;
  bsum =
      rt_powd_snf(10.0, (alpha_idx_0 * log10(n) + alpha_idx_1 * log10(L * B)) +
                            alpha_idx_2) *
      100.0;
  /*  判断估计质量是否偏差过大 */
  if (fabs(bsum - original_mass) / bsum > 0.5) {
    original_mass = bsum;
    original_Modal_Mass1 = bsum * C / n;
  }
  *eachmass = original_mass / n;
  i = y_val2->size[0] * y_val2->size[1];
  y_val2->size[0] = 1;
  y_val2->size[1] = (int)n;
  emxEnsureCapacity_real_T(y_val2, i);
  y_val2_data = y_val2->data;
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val2_data[i] = 0.0;
  }
  if ((int)n - 1 >= 0) {
    for (i = 0; i < 101; i++) {
      y[i] = a[i] * n;
      b_y[i] = phi[3 * i + 1];
    }
  }
  if (loop_ub_tmp - 1 >= 0) {
    spline(y, b_y, c_y, t1_coefs);
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val2_data[i] = ppval(c_y, t1_coefs, (double)i + 1.0);
  }
  i = y_val3->size[0] * y_val3->size[1];
  y_val3->size[0] = 1;
  y_val3->size[1] = (int)n;
  emxEnsureCapacity_real_T(y_val3, i);
  y_val3_data = y_val3->data;
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val3_data[i] = 0.0;
  }
  if ((int)n - 1 >= 0) {
    for (i = 0; i < 101; i++) {
      y[i] = a[i] * n;
      b_y[i] = phi[3 * i + 2];
    }
  }
  if (loop_ub_tmp - 1 >= 0) {
    spline(y, b_y, c_y, t1_coefs);
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    y_val3_data[i] = ppval(c_y, t1_coefs, (double)i + 1.0);
  }
  /*  计算刚度矩阵和 Modal Mass */
  i = b_I->size[0] * b_I->size[1];
  b_I->size[0] = (int)n;
  b_I->size[1] = (int)n;
  emxEnsureCapacity_real_T(b_I, i);
  I_data = b_I->data;
  for (i = 0; i < b_loop_ub_tmp; i++) {
    I_data[i] = 0.0;
  }
  for (i = 0; i < loop_ub_tmp; i++) {
    I_data[i + b_I->size[0] * i] = 2.0;
  }
  i = c_d->size[0] * c_d->size[1];
  c_d->size[0] = (int)(n - 1.0) + 1;
  c_d->size[1] = (int)(n - 1.0) + 1;
  emxEnsureCapacity_int8_T(c_d, i);
  d_data = c_d->data;
  for (i = 0; i < c_loop_ub_tmp; i++) {
    d_data[i] = 0;
  }
  for (i = 0; i < nv_tmp_tmp_tmp; i++) {
    d_data[(i + c_d->size[0] * i) + 1] = -1;
  }
  i = d_d->size[0] * d_d->size[1];
  d_d->size[0] = (int)(n - 1.0) + 1;
  d_d->size[1] = (int)(n - 1.0) + 1;
  emxEnsureCapacity_int8_T(d_d, i);
  b_d_data = d_d->data;
  for (i = 0; i < c_loop_ub_tmp; i++) {
    b_d_data[i] = 0;
  }
  for (i = 0; i < nv_tmp_tmp_tmp; i++) {
    b_d_data[i + d_d->size[0] * (i + 1)] = -1;
  }
  m0 = sparse(n, n, n, n, t2_d, t2_colidx, t2_rowidx, &acol_size);
  if (b) {
    i = c_I->size[0] * c_I->size[1];
    c_I->size[0] = b_I->size[0];
    c_I->size[1] = b_I->size[1];
    emxEnsureCapacity_real_T(c_I, i);
    C_data = c_I->data;
    for (i = 0; i < b_loop_ub_tmp; i++) {
      C_data[i] = (I_data[i] + (double)d_data[i]) + (double)b_d_data[i];
    }
    sparse_plus(c_I, t2_d, t2_colidx, t2_rowidx, m0, acol_size, b_I);
    I_data = b_I->data;
  } else {
    binary_expand_op(b_I, c_d, d_d, t2_d, t2_colidx, t2_rowidx, m0, acol_size);
    I_data = b_I->data;
  }
  emxFree_real_T(&c_I);
  emxFree_int32_T(&t2_rowidx);
  emxFree_int32_T(&t2_colidx);
  emxFree_real_T(&t2_d);
  emxFree_int8_T(&d_d);
  emxFree_int8_T(&c_d);
  acol_size = y_val->size[1];
  c_i = b_I->size[1];
  i = b_C->size[0] * b_C->size[1];
  b_C->size[0] = 1;
  b_C->size[1] = b_I->size[1];
  emxEnsureCapacity_real_T(b_C, i);
  C_data = b_C->data;
  for (i = 0; i < c_i; i++) {
    m0 = i * b_I->size[0];
    C_data[i] = 0.0;
    for (b_k = 0; b_k < acol_size; b_k++) {
      C_data[i] += y_val_data[b_k] * I_data[m0 + b_k];
    }
  }
  emxFree_real_T(&b_I);
  emxFree_real_T(&y_val);
  C = 0.0;
  m0 = b_C->size[1];
  for (i = 0; i < m0; i++) {
    C += C_data[i] * c_data[i];
  }
  emxFree_real_T(&b_C);
  emxFree_real_T(&c);
  *k = Ki * original_Modal_Mass1 / C;
  /*  返回 Modal Mass 和矩阵 */
  bsum = 0.0;
  m0 = y_val2->size[1];
  for (i = 0; i < m0; i++) {
    b_d = y_val2_data[i];
    bsum += b_d * *eachmass * b_d;
  }
  emxFree_real_T(&y_val2);
  b_expl_temp = 0.0;
  m0 = y_val3->size[1];
  for (i = 0; i < m0; i++) {
    b_d = y_val3_data[i];
    b_expl_temp += b_d * *eachmass * b_d;
  }
  emxFree_real_T(&y_val3);
  Modalmass[0] = original_Modal_Mass1;
  Modalmass[1] = bsum;
  Modalmass[2] = b_expl_temp;
}

/* End of code generation (Modalmass_estimation_newchangefinal1219.c) */
